Mode coupling theory for multi-point and multi-time correlation functions 
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We present a theoretical framework for higher-order correlation functions involving multiple times 
and multiple points in a classical, many-body system. Such higher-order correlation functions have 
attracted much interest recently in the context of various forms of multi-dimensional spectroscopy, 
and have found an intriguing application as proposed measures of dynamical heterogeneities in 
structural glasses. The theoretical formalism is based upon projection operator techniques which 
are used to isolate the slow time evolution of dynamical variables by expanding the slowly-evolving 
component of arbitrary variables in an infinite, "multi-linear" basis composed of the products of slow 
variables of the system. Using the formalism, a formally exact mode coupling theory is derived for 
multi-point and multi-time correlation functions. The resulting expressions for higher-order correla- 
tion functions are made tractable by applying a rigorous perturbation scheme, called the A'^-ordering 
method, which is exact for systems with finite correlation lengths in the thermodynamic limit. The 
theory is contrasted with standard mode coupling theories in which the noise or fiuctuating force 
appearing in the generalized Langevin equation is assumed to be Gaussian, and it is demonstrated 
that the non-Gaussian nature of the fiuctuating forces leads to important contributions to higher- 
order correlation functions. Finally, the higher-order correlation functions are evaluated analytically 
for an ideal gas system for which it is shown that the mode coupling theory is exact. 



PACS numbers; 61.20.Lc, 05.20.Jj, 05.40.-a. 
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I. INTRODUCTION 

In a glassy system where structural frustration pre- 
vents relaxation to equilibrium, dynamical properties of- 
ten demonstrate complicated time dependence Q. For 
instance, in dense colloidal systems, at a given time some 
regions of the complex fluid are essentially static and 
crystalline, while the dynamics in other regions exhibit 
behavior which is associated with fluids. In these sys- 
tems, structural rearrangement occurs through relatively 
rapid, collective, string- like motions Furthermore, 
at later times, a region of the fluid which previously ap- 
peared crystalline may exhibit fluid-like properties. Such 
heterogeneous behavior is characteristic for structural 
glasses and super-cooled complex fluids To de- 

scribe this behavior, it is natural to examine how the lo- 
cal densityof the liquid is correlated over various spatial 
domains |]lO|,0, or, when one is more interested in the 
different time scales of slow global changes of structure 
and the local decay of correlations, multiple time corre- 
lation functions |l2| . Both types of higher-order correla- 
tion functions have recently been proposed and used as 
measures for characterizing "dynamical heterogeneities" 
in structural glasses. 

A number of experimental probes of examining de- 
tailed dynamical features taking place on various length 
and time scales in glassy systems have emerged over the 
last few years. These new approaches have the poten- 
tial to provide extremely useful information on how col- 
lective motions of the system are correlated to specific 
statistical features of the dynamics such as the distribu- 
tion of time scales of fluctuations, the length scale and 
size-distribution of solid-like clusters, and cage structural 



relaxation rates. One approach to probe the nature of dy- 
namical heterogeneities is based on single molecule spec- 
troscopy techniques |13-16| in which the environment of 
a one molecule is tracked over a period of time. The 
technique allows the information of not only the distribu- 
tion of heterogeneous environments but also the explicit 
reorganization times which are present in the system, 
since the individual measurements are not statistically 
averaged. A somewhat different experimental approach 
is based on multi-dimensional NMR | |r^ , p^ and non- 
resonant non- linear Raman spectroscopy |19|-pl[. The 
response function in these experiments can be related to 
higher-order correlation functions using response theory 
||,||. 

Given the interest in multi-time and multi-point cor- 
relation functions, the need for a theory that accurately 
predicts these quantities is clear. Surprisingly, there has 
been relatively little work along these lines and the lit- 
erature is not nearly as extensive as it is for ordinary, 
two-time, two-point correlation functions such as the dy- 
namic structure factor. Although there have been sev- 
eral recent microscopic theories for the off-resonant fifth- 
order response function for simple liquids |p^ , p5| , little 
work has been done on constructing microscopic theories 
for general higher-order correlation function since the ki- 
netic theories of De Schepper and coworkers who 
attempted to extract the non-analytic density contribu- 
tions to the Burnett coefficients in hard sphere liquids. 

A common approach to describing the dynamical prop- 
erties of liquids and complex fluids at long times is based 
on the generalized Langevin equation pTf . The basic util- 
ity of the generalized Langevin equation depends on the 
assumption that the long-time behavior of an arbitrary 
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dynamical variable of the system can be written in terms 
of the dynamics of a specific set of slow modes. This slow 
behavior can be isolated by extracting the projection of 
the dynamical variable onto the slow modes, which effec- 
tively form a basis set for the long time behavior of the 
system. This approach has been successfully applied to 
describe relaxation and simple time correlation functions 
in a wide variety of condensed phase systems. In the 
context of simple liquid systems, it was initially assumed 
that only the linear densities of conserved variables of 
the system composed the set of slow variables of the sys- 
tem p^-|30|. However, it was discovered that this basis 
set was insufficient |3l|] to describe the complex, non- 
exponential relaxation of simple time correlation func- 
tions observed in molecular dynamics simulations [^ . 
The subsequent observation that theories incorporating 
multi-linear products of the linear densities in the basis 
were capable of yielding the correct asymptotic long-time 
behavior of simple correlations [ p3[|34| led to the devel- 
opment of "mode-coupling" theories. The seminal work 
of Kawasaki | |35| ] , who proposed that the linear Langevin 
equation be replaced by a non-linear version in which the 
fluctuating forces obey Gaussian statistics, sparked the 
later development of kinetic mode-coupling theory mod- 
els of dense liquids |3^-|4C|]. At roughly the same time, 
Ronis used the framework of Kawasaki model cou- 
pling theory to formulate a theory of higher-order corre- 
lation functions in which the multi-linear slow variables 
forming the basis set for long-time evolution in the sys- 
tem obey Gaussian statistics. Although Ronis' theory 
contains a number of assumptions, it is the assumption 
of Gaussian statistics which leads to clear inconsistencies 
in the predictions for higher-order correlation functions. 
Unfortunately, the Gaussian assumption is fundamental 
in the Kawasaki formulation of mode coupling theory, 
and is difficult to generalize. 

The purpose of this paper is to provide a solid the- 
oretical framework to calculate multi-point and multi- 
time correlation functions. It is based on the mode cou- 
pling theory obtained by a projection operator formal- 
ism as developed by Oppenheim and coworkers 
p^-{47[. By careful consideration of how to consistently 
identify fast and slow behavior in time correlation func- 
tions, we derive expressions for multi-point and multi- 
time correlation functions in terms of simple time corre- 
lation functions. Since the basis set for the slow modes 
is infinite, an infinite number of terms arise in the ex- 
pressions for the higher order correlation functions. It 
is demonstrated that if the system has a finite correla- 
tion length, the infinite series can be truncated by ap- 
plying the so-called A'^-ordering perturbation expansion 
method, which is exact in the thermodynamic limit. The 
use of this perturbation scheme circumvents the need 
to assume that the basis set obeys Gaussian statistics. 
Based on this method, the leading order expressions and 
first mode-coupling corrections for higher order correla- 
tion functions are presented. Finally, it is shown how the 
formalism applies for an ideal gas system for which it is 



demonstrated that the theory yields the exact result for 
simple, multi-point, and multi-time correlation functions. 

II. TWO-TIME CORRELATIONS 

A. The system and slow variables 

Consider a classical system composed of A^ point par- 
ticles in which the momentum and position of parti- 
cle i are denoted by Pi and r.^, respectively. Given 
the Hamiltonian Ti, a function B{T) of the phase point 
r = (ri, • • • , rjv, pi, • • • , PAf) evolves according to 

Bit) ^ {Bit),H} = CBit), 

where {, } denotes the Poisson bracket and £ is the Li- 
ouville operator for the system. The evolution equation 
can be solved formally as B{t) = exp{£t)B{0), where 
here and below B{t) is taken to denote B(r{t)). 

Typically, each dynamical variable of the system can be 
separated into slowly- varying and quickly- varying parts. 
We will assume that the time-dependent correlation func- 
tions of the quickly- varying components decay to zero on 
a short (called "microscopic") time scale Tm, whereas cor- 
relation functions of the slowly varying part decay on a 
longer time scale Th- Hence, at long times t ^ Tm, the de- 
cay of an arbitrary correlation function can be described 
by the decay of its slow component. In what follows, we 
postulate that the slowly varying part of an arbitrary dy- 
namical variable is an analytic function of a set of slow 
variables of the system. In this sense, the slow variables 
form a basis set in which to represent the long-time be- 
havior of the system. 

To identify slow variables in such a system, one con- 
siders all the conserved quantities, which can be taken 
together in one column vector A, with components A"'. 
When these quantities are extensive, they can be ex- 
pressed as a sum over contributions from the individual 
particles, 

N 

yi(r) = ^a,(r), (1) 

which leads to the identification of the densities as the 
local version of A, 

A(r;t) = ^a,(<)<5(r-r,(t)), (2) 

3 

with Fourier components 

A(<) = E«.(*)^"'""^'^*^- (3) 
j 

For the case of a simple fluid of N point particles, the 
extensive slow variables are the number density, momen- 
tum density, angular momentum density and the energy 
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density, and in Eq. (g), one would use the microscopic 
expressions for those quantities. For point particles, the 
angular momentum density can be expressed in terms of 
the momentum density and need not be included in A. 

The small wave- vector components of the densities 
correspond to large length scale fluctuations, and these 
are expected to evolve slowly since their time derivatives 
are proportional to |k|. On the other hand, large wave- 
vectors with magnitudes beyond some cut-off value kc 
correspond to small length-scale fluctuations and have 
large time derivatives provided the system is not too 
dense. Thus, one can identify the slow variables to be 
composed of the Fourier-transform of densities of con- 
served variables whose wave- vector arguments are re- 
stricted have magnitudes less than fc^. 



B. Generalized Langevin equation 



The distinction between fast and slow behavior can be 
made more precise on the level of correlation functions, 
following the projection operator method as used by Op- 
penheim and others p^-^. 

To isolate the slow part of a correlation function, one 
assumes it to be an analytic function of the slow vari- 
ables The non-linear dependence of the time cor- 
relation functions on the slow variables can be incorpo- 
rated within the framework of projection techniques by 
employing a basis of non- linear functions of Ak, the so- 
called multi-linear basis. 

Using the ensemble average (• • • ) as an inner product, 
we define a projection operator Vi that projects onto the 
deviations ^k = — (^k) of the slow variables from 
their equilibrium value, as 



Qo 



— ^k 

_ /ja 



A'' 



Aa 
^k 



iro + ri)At 



A'' 
q^q 



(5) 



Qai,a2,.--0'Ti 
ki-k',k2,... ,k„ 



J=0 



/.ai Aa2 Aa^ 

"^ki-k'^k2 "^k„' 



where k' — X]r=2 ^'^'^ projection operators Vj are 
defined as 



E 



{XQs 



aa 



(6) 



In this notation, a Greek index denotes a set of wave- 
vectors and hydrodynamic indices, and |a| denotes the 
number, or mode order, of hydrodynamic indices in a set 
a. For simplicity of notation, the full notation will often 
be written Qi, Q2, Q3, etc., for Qa when |q;| — 1, 2, 3, 
.... Furthermore, non-hatted and hatted Greek indices 
will always have the same mode order, i.e. \a\ — \a\, 
but represent different sets of wave- vectors and/or hy- 
drodynamic indices. Finally, in Eq. (^) Kaa — {Q&Qa), 
and the -product now includes a sum over all hydro- 
dynamic indices and wave-vectors in the summation in- 
dices (a and d). To make sure the contributions from the 
same component in a sum over a are not over-counted, 
one has to divide by the number of ways the indices can 
be rearranged in a. 

The projection terms in Eq. (||) force the set Qa to 
be orthogonal in mode order so that {QaQp) = unless 
|a| = \f3\. This property will be not only convenient but 
very important in the subsequent analysis. By assump- 
tion, the long-time behavior of an arbitrary variable C 
can be isolated by the projection operator, 



where 



ViX = (XAk) * i^k ' * ^k, 



^k = (^kij;) 



(4) 



denotes the normalization. The "*" -product in Eq. (^) 
includes a sum over components of the column vector 
Ak - the hydrodynamic indices - as well as over wave- 
vectors k. For a translationally invariant system, all 
wave- vectors in an average must add to zero, so only 
when k equals the sum of wave-vectors in X will there 
be a contribution in Eq. (^). 

Using only the linear projection Vi will not give a 
proper mode coupling theory since the long-time depen- 
dence of dynamical variables due to multi-linear orders of 
the densities is not extracted by the projection operator 
p6|-^. It is therefore convenient to define an orthogo- 
nal, multi-linear basis as 



3=0 



(7) 



where we have used the convention that repeated Greek 
indices imply a -product and a summation over mode- 
order. This notation will be used throughout this article 
unless stated otherwise. 
Applying operator identity 

Jo 

to the evolution equation, one easily obtains. 



f e^'^'-^^VCe^^^P^dT, 
Jo 



(8) 



where V± — 1—V and C±_ = V^C. We apply this operator 
to Qa, and write out V to get 
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•J 

Mo,fi{T)Qp(t-T)dT + dp^{t), (9) 



which is the generahzed Langevin equation, where 

M^0{t) ^ [2<5(r)(Qagp - {Mr)<Pl)]Kjl, (10) 

and the fluctuating force (f>a{t) is defined by 

0„(t) = e^^*£^g„ = e^^*7'_iQ„. (11) 

Under the assumption that V projects out all the slow 
behavior, the matrix Mcfj{T) is "fast" in the time vari- 
able T in the sense that, for a microscopic time scale Tm 
much smaller then the hydrodynamic time scale Th of the 
slow variables, 

{Mr)B)^ {Mt)B)dtS{T)+0{Tra/Th). (12) 

JO 

It then follows that the main contribution in the integral 
in Eq. comes from small r, and we can approximate 
Eq. (^ by one which is local in time, 

Qa{t) ^ MapQp{t)+4'c.{t), (13) 

where the instantaneous matrix is given by 



(14) 



It will become clear in section |IIE| that the assump- 
tion of a fast decaying memory kernel in Eq. (^) is not in 
contradiction with the slow memory kernel that occurs 
in linear mode coupling equations like those used in the 
study of the glass transition Q . 

Defining the time correlation functions of the basis set 
to be 



M„^= / M^p{t)dt. 



(15) 



and using Eq. (^, one obtains the simple expression for 
the time correlation function, 



Ga/3(i) = / M^s{T)Gsp{t-T)dT, 



(16) 



where we have used the fact that {4'a{t)Q*p) = by con- 
struction. 



In the instantaneous approximation [corresponding to 
Eq. (|l])], Eq. (|l|) yields G^p{t) = M^sGsp{t), which 
can be integrated to obtain 



G a/lit) 



[e Ja/3- 



(17) 



If we do not wish to make the instantaneous approxima- 
tion, it is easiest to look at the Laplace transform of the 
time correlation function 



G^p{z)~ / G^p{t)e-'' dt. 
Jo 

Since the time convolution in Eq. (|l6| ) is a simple product 
in Laplace space, we have 

Gaf^iz) = {Qa{z)Ql)K7^, = [zl-M{z)]~' (18) 



where the a (3 element of the inverse is meant, and 



M^piz) = [(Q„gt 



haiz) 



0" P0 

Note that one is frequently interested in just the part of 
Gap{t) that involves the linear variables, corresponding 
to such readily observable physical quantities as the time 
correlation functions of the linear density, momentum 
and energy. These are given by Gii(i), the |a| = |/?| = 1 
sub-block of the infinite-dimensional matrix Gapit). 

In the hydrodynamic limit in which the magnitude of 
the wave-vectors k is small, the wave-vector can be used 
to perturbatively order various contributions to Map{z). 
Noting that each time derivative in Eg . (p^ brings down 
a factor of k, the first term in Eq. ([lO|), called the En- 
ter term, is of order k, whereas the second, or dissipative 
term, is 0(|kp). On the basis of these arguments, one 
might think that the dissipative term in Mq^(z) can be 
neglected in the hydrodynamic limit. However, exam- 
ining the form of the densities, it is easy to see that in 
fact the Euler term is imaginary and can only give rise to 
oscillatory behavior^. It is therefore clear that the dissi- 
pative term must be included in order to obtain solutions 
to Eq. (O) which are well-behaved in the long-time limit. 



C. TV-ordering 

Since two-time correlation functions at a particular 
mode-order are given by specific blocks of the inverse 
of an infinite dimensional matrix, even the calculation 
of simple time correlation functions requires an infinite 
number of multi-linear modes. Thus, the results of the 
previous section are formally exact, but not very useful. 

The evaluation of the time correlation functions is 
greatly facilitated by applying cumulant expansion or A^- 
ordcring techniques. The A^-ordering scheme was first 



^The Euler term can give rise to Gaussian relaxation if one 
has an infinite number of slow modes, see Sec. ^ 
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introduced by Machta and Oppenheim [Q as an exten- 
sion of Van Kampen's inverse system size expansion {fl 
expansion) and developed further in Ref. |Q. It is 
essential in obtaining the correct Stokes-Einstein law for 
a Brownian particle using mode coupling techniques p7| . 

In the iV-ordering approach, one assigns to each cumu- 
lant of an average appearing in the equations an order of 
N (the number of particles). The starting point is to 
consider cumulants, which we denote by "((• • • ))". For a 
product of linear densities, the cumulant expansion con- 
sists of all possible ways of combining the densities into 
groups, i.e., (A) = {{A)), (AB) = {{AB)) + ((A)) ((i?)), and 
so on. The assignment of TV-orders is based on the obser- 
vation that each cumulant containing n linear densities 
is of order iV(^/a)'^'-"~^\ where a is the average distance 
between particles ^ij. Hence, the requirements for the 
expansion method to be meaningful are that the system 
should have a finite [0{N'^)] correlation length ^, and 
that the integrated densities should be extensive. For an 
extension to non-extensive quantities like tagged particle 
densities, see Ref. |^ . 

For instance, the cumulant of a linear-linear corre- 
lation function of the number density A^k, defined as 



N 



, can be formally written as 



((iVkA^k)) = (^kiVk) - (A^k)(A^k>4o 



[{N) + (N{N~l}e 



ik-(ri-r2; 



>](i-M 



kO- 



For k = 0, the cumulant expansion of {N^N^) is given by 
{{N — (N))'^), which, in the grand canonical ensemble, is 
of order (N). For k ^ 0, the expression is proportional to 
N as well, because particles beyond the correlation length 
will not contribute to the average of exp(ik • (ri — r2)), 
so that 



{N{N - 1) exp(ik • (ri - ra))) cx N{N - 1) 



V 



0{N). 



The iV-ordering of various quantities has been dis- 
cussed extensively elsewhere |^,^,^, so here we just 
state a few properties of the procedure that enable us to 
estimate A^-orders of relevant quantities. 

One important basic property of averages involving Qa 
is that when one has an average of the form {Qa{t)BQl), 
where B can be a product of Qi again, it can be factored 
into a term of order iV^ , 



{Qait)BQZ*) « {Q^^i{t)B){QlUt)Ql*)S: 



k,k 



(19) 



plus terms involving lower powers of N. That is, the lead- 
ing A^-order term in the cumulant expansion can be found 
by equating any wave- vector kj from the set a with k. In 
Eq. (|l^), a — 1 denotes the set a with kj and aj removed. 
Hence, means that the expression is correct up to 
higher orders of A^ and only for the case where appropri- 
ate wave- vectors have been equated. We can reduce such 
an expression further by equating wave- vectors in B with 



those in the set a — 1. An important point, established 
in Ref. |^] , is that equating wave- vectors within the set 
a does not increase the A^ order due to the subtraction 
terms in their definition Eq. (||). Thus the orthogonal- 
ization procedure used to construct the multi-linear basis 
plays an important role in the proper A'^-ordering of cor- 
relation functions of multi-linear densities. 

Using the ordering scheme, one obtains Kaa = 0{N^), 
where j is the number of matched sets of wave- vectors in 
a and d, and hence the A^-ordering of its inverse is 



k: 



0{N- 



(20) 



irrespective of matching of wave- vectors. In Ref. |4J], it 
was shown that 



Oil) 

0(iVl"l- 



-W) 



if |a| > 
if |a| < 



with the same A^-ordering holding for Map. A similar 
result can be obtained for Gafj{t) 



0(1) 
0(7Vl' 



if |a| > 
if |a| < 



(21) 



D. Multi-point correlations 

We define a multi-point correlation function to be any 
Gapit) for which \a\ or \(3\ is greater than 1. The multi- 
point correlations functions therefore contain several fac- 
tors of the linear densities Ak, and involve more than 
one wave-vector. Hence, in position space, the multi- 
point correlation functions involve at least three spatial, 
or two relative, position coordinates, justifying the name 
multi-point correlation functions. 

As Map{z) is 0(1) if a and /3 have all wave- vectors 
matched, and of lower order otherwise, the leading A^- 
order contributions to Gap{t) can be obtained by ex- 
panding Eq. (18) in the wave- vector-diagonal compo- 
nent of Map{z), denoted by = Maa'{z)5a' p- Here 
primed Greek indices will always have the same set of 
wave- vectors as as their unprimed variant, but not nec- 
essarily the same hydrodynamic indices. Defining the 
off-diagonal part as M°p{z) = Mapiz) — Maa'{z)Sa'p, 
we get from Eq. ( |l8| ) 

Gapiz) = [zt- M'^iz) ^ M°iz)]-l 

^ [{z]l-M'^(z)}{]l-G(z)Ar(z)}];^ 

= [l-Giz)M"iz)]-l,Gp,piz), (22) 

where we have defined the diagonal in wave-vector 

Gao.'iz)^[zt-M^{z)]-l- (23) 
We can therefore expand Eq. (p2h as 
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+Ga,^,{z)M°,^{z)G^y{z)M°,i^,{z)Gfj,p{z) 
+ ... (24) 

Note that when a term is of lower order in N, this 
does not mean it can be neglected in the thermodynamic 
limit. The leading iV-order is found from equating wave- 
vectors, whereas the next order in N often comes with an 
un-restricted summation over wave- vectors. If M denotes 
the number of (slow) wave-vectors in the system, these 
correction terms are 0{M/N). The number of wave- 
vectors M grows with the system size. Therefore such 
mode coupling corrections survive in the thermodynam- 
ics limit. But when the mode coupling parameter M/N is 
small, as is typically the case for simple liquids at moder- 
ate densities, it can be treated as perturbation parameter 



where by definition, Ggp' has the same form as the right- 
hand side of Eq. (p^, with a replaced by 6 and the re- 
striction that none of the wave- vector sets in the expres- 
sion are equal to a. Being of that form, we can repeat 
this procedure for Gg'^ to obtain. 



Gaf3 — Gaa'^a 'P + Gaa'M°gG 



(a) 
■5/3 



When the procedure is continued ad-infinitum, we find, 

Gap = Gaa'Sa'13 + Gqq' Af °/^/ G^,^ 
+Gaa' M°gG'gJ Mg,0, Gp,p ^ 



In Eq. ( p4| ) , note the appearance of Gaa' , defined in 
Eq. (|2^). Although this is diagonal in wave- vector, it is 
not the diagonal part of Gap , as the terms following the 
first term in Eq. (|^ give contributions for a and /3 diag- 
onal. As was shown in Ref. p6|] , in the thermodynamic 
limit, the diagonal part of Gap, denoted by Gaa', can be 
factored as 



Gaa' (t) ~ 



k,k' 



(25) 



where aj and a'^^ are hydrodynamic indices from a and 
a', and and k^^. the respective wave- vectors. The 
summation is over all permutations a of the indices in 
a'. This factorization is obtained also by cumulant ex- 
pansion under the assumption that there is a finite time 
dependent correlation length. 

Thus, if we were able to express Eq. ( p4[ ) in terms 
of Gaa' instead of Gaa', we could combine that with 
Eq. ( p5| ) to get an expression for any multi-point func- 
tion in terms of in the two-point correlation functions Gn 
and vertices M°^. Eq. ( p4[ ) can be re-expressed in this 
desired form by the following re-summation of terms: We 
write Gap — Gaa'Sa'p + G°^p, and use the Dyson form of 
Eq. (H): 

Gap = Gaa'^a'P + Ga-fM°p,Gp'p 

for the off-diagonal part to get 

Gap = Gaa'^a'P + Ga^M°p,G^^,^^, 

where the superscript (a) means that /3 is restricted to 
not have the same wave-vector set as a. Iterating this 
equation yields 

Gap — Gaa'Sa'P + G aa' M^t f^i G^p}^ 

+Gaa' M°,^,G[^,}^M°i^,G^p,^p + • ■ • 

= Gaa'^a'P + Gaa' M°gG^i^ , 



This expression resembles that for Gaa' in Eq. (p^) , with 



P = a' and G^^' replaced by G, 



'77 
'.(a,...) 



Furthermore, the 



definition of G^y " also resembles the definition of G^^' , 
but now with restrictions on all wave- vector sets. In fact 
the wave-vector restrictions can be relaxed in the ther- 
modynamic limit since the restrictions remove only one 
term out of the sum over intermediate wave-vector. Rel- 
atively speaking, the difference between the series with 
restricted and unrestricted sums is of order 0{1/N), so 
the restriction on the intermediate wave-vectors is negli- 
gible in the thermodynamic limit, and we can write 



G 



(a,...) 
PP' 



G 



(a,. 

PP' 



where G^^^i"^ is the full correlation function which is 
diagonal in wave-vector and in which the set of wave- 
vectors in P differs from those in sets in a, ... . 

Thus Eq. (p^), which has the non-physical G, can be 
replaced by the expansion 



Gap — Gaa' 5a' P 



- Gaa' M°, f^iG^pijj 



-Gaa' M°gG''gJ ^S'P' G 



ia,S) 
P'P 



-tUaa'MaS^SS' ^^^S'l^-y^' ^^^-f'p'^p'p 



(26) 



which involves the full correlation function Gaa' which is 
diagonal in wave-vector. In this expression none of the 
intermediate wave-vector sets are allowed to be equal. 
Using Eq. (p6|), we can write to leading iV-order 

Gi2{z) = Gii(z) * Mi2(z) * G22{z), (27) 

where, from Eq. (^), G22{t) is given by, 

^k-q',q';k-q,q(0 ~ G'k-q(^) * G]^ {t) (5q,k-q' 

+ Gll_q(t)oGll(0<5qq', (28) 
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where superscripts like 22 and 11 are a reminder of the 
mode orders of the arguments. In order to faciUtate writ- 
ing tensor products, we have introduced the following no- 

tational symbols for products of tensors of rank 2 {A"-'''), 
rank 3 (yl"-^^ if b and c belong to the same set a or A"-''''' 
if a and b belong to the same set), or rank 4 (such as 

{A ■ Bf" = A^'^B^'" 
[A ■ B)"''^ = A"'''^'^ B'^'f''^ 

{Ao BY'^'-"''^ = A'^'-'^B''''^ 

^j^,^yi,b;c,d^^a;dBb;c^ (29) 

where repeated labels are summed over. 

The two terms in the expression above for G22 turn 
out to yield the same contribution to G12 in Eq. (27). 
However, this summation has a pre- factor of 1/2 from 
the number of ways the indices can be interchanged in 
the product. Using Eq. {M) in Eq. (|27|), we obtain 



'^kik-q.qlO — / (i - r) • Afk.i^_q q 

Jo 

: [Gil^{T)oGl\r)]dT, (30) 

where the time convolution arises from the inverse 
Laplace transform. 

At this point, the necessity of including multi- linear 
modes is readily apparent. For example, using the defi- 
nition of the multi- linear basis set, we can write 

(QimiQl) - (Qi(t)QD */^ri' * (QiQlQl) 

+Gi2{t)*K22, (31) 

and note that the second term would have been absent 
if the bi-linear modes Q2 in the basis set hadn't been 
included even though the A^-ordering of this term is the 
same as the first term. Another interesting point is that 
the first term is not present if the subtractions in the 
definition of Q2 in the basis set are not included. If one 
assumes that the Qi are Gaussian random variables, then 
correlations of the form {QiQlQD vanish. However, in 
dense fluids the linear densities clearly do not obey Gaus- 
sian statistics since static correlations such as {QiQlQl) 
involve configurational averages over the triplet distribu- 
tion function and are not negligible. Note that the second 
term involves a time convolution in Eq. ( pO[ ) and can be 
expected to have quite different behavior from the first 
term, which is proportional to an ordinary time correla- 
tion function. 



E. Renormalization of the propagator 

In this section we focus on the linear correlation func- 
tion Gil itself. Eq. (p3) for Gii(z) reads: 



Gil — Gil + GiiM°^Gaa' M°,^Gii 

+G11 Gaa' M°, 13, G 010 Mpj^Gii 

+GiiM°^Gaa' M°, p,G0' 0MI^^Gyj' M°,j^Gii 
+ ... 

where for brevity we omitted the z argument. In the 
summations over a, one can isolate all the terms with 
mode order 1 to obtain, 

Gil — Gil + Gii(M°^Gc(C(/M°,]^)Gii 

+Gii(M{'^Gq,q'M°,^,G/3'^M^j)Gii 
+G11 (^Adi^Gaa' M°,0,G0i0MI^^G^^' Myi'^Gii 
+ G 1 1 ( M{'q, G Q, Q ' M ° , 2 ) G 1 1 ( G^^i M°, 1 ) G 1 1 
+ ... , 

where the summation over repeated indices here start at 
mode order 2. Rearranging the terms, this can be written 
as 

Gil = Gil + Gil * Sii * Gil 

-fGii * Sii * Gil * ^11 * Gil + • . • , (32) 

with 



(33) 



Eq. ( ^ ) can be re-summed as 

Gii(z) = [z]l-Afii-Eii(z)]-i (34) 

where the inverse is taken on the 11 sub-block level. 
Eq. (Q) can be utilized to extract the complicated long- 
time dependence that arises in the memory functions for 
the generalized Langevin when only linear densities are 
included in the projection V . In the theory of liquids, the 
Laplace transform of these memory functions are gener- 
alized transport coefficients which reduce to the Green- 
Kubo expressions in the limit of small z and k. Eq. (jsj) 
can be cast in the form in which the full generalized trans- 
port coefficients are expressed as a sum of bare transport 
coefficients and the I]ii(z) terms, which renormalize the 
bare coefficients and account for the complicated mem- 
ory effects observed in dense liquids |^^. In the limit 
in which the energy density is neglected in the basis set 
and only bilinear modes are included in the multi-linear 
basis set, the idealized and extended mode coupling the- 
ory models |51| of the glass transition can be obtained 
P6|,PI. 
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Finally, we mention that in Ref. j4q| it was shown that 



by rearranging the terms, along the lines of Sec. [I E , Gaa' 
in the Qi can be replaced by the real diagonal in wave- 
vector Gaa', with the restriction in summations over a, 
/3 etc, that none of their the wave- vector set are identical. 
This diagonal Gaa' will factor again as in Eq. (p5|). Then, 
Eq. dH), Eq. @ and Eq. (H) lead to a self-consistent 
equation for Gn. 



III. MULTI-TIME CORRELATIONS 

A. Separating slow and fast behavior in the context 
of multiple times 

We now turn to multiple time correlation functions like 
{Q{t2 + ti)Q{ti)Q)- At first glance, it appears attractive 
to use the formal formal solution of Eq. (0) , 



Msne - <! - / / {4>5{T2)Qi3(t)l{-ri))dTidT2 

JQ Jo 



Q<.{t)^Gasit)Qs{0)+ Ga5{t-T)Mr)dT, (35) 
Jo 

to get an exphcit expression for (Qa{ti + t2)Q p{ti)Q^) ■ 
After all, inserting Eq. (^5|) yields 

{Qo.{ti^t2)Qp{ti)Q*) 



tl+t2 



dT2 



dri Ga-niti +t2 - T2)Gf3s{ti - n) 



{QnQ5Q;HS{T2)S{Ti) + {MT2)Mri)Q; 



(36) 



It is tempting to assume that the expression in Eq. ( pq ) 
involving the fluctuating force behaves as S(t-j)S(t-\ ) at 
long times in an analogous fashion to Eq. (|lj), leading 
directly to a local equation in time, 

{Qo.itl+t2)Ql3itl)Q;) = Gariitl+t2)Gpsiti)MrjS-y, 

(37) 

where M^^^ is related to the infinite time integral of 
[{QsQsQ;)^SiT2)SiTi) + (,^^(t2)05(ti)O;)]. However, 
things are not this simple as can be seen by using time 
translation invariance to write {Qa{ti + t2)Qg{ti)Q%) = 
{Qa{t2)Qi3Q*^{-ti))- Applying again Eq. (||), we obtain 

{Qc.{t2)QpQ^-h)) 

dT2 I dTi Garj{tl + h - T2)G*s{~ti - Ti) 

Jo 

{Qr,Qf3Q*s)'^S{T2)S{T,) + {c^,{T2)Qp(l>Kri)) , 

of which the local time version would be 

7 _ 

{Qa{t2)QpQ*^{ — ti))Kxi^ = Garj{t2)M^I3sGsj{ti) (38) 

where we have used the time-translation property, 
^api^)^pp = Gpa(-t)K&a (scc Appendix]^), and de- 
fined 



(39) 



Clearly, Eq. ( p^ ) and Eq. ( pTj ) are in contradiction since, 
in general, they will have different time behavior. 

The question of which (if either) instantaneous form 
is approximately correct can be resolved by noting that 
products of the fluctuating force (/)(t) cannot always 
be treated as "fast", nor as Gaussian random vari- 
ables. This observation has been noted previously by 
Schramm and Oppenheim p5| who considered the quan- 
tity ((/)(Ti)(/)(T2)(/)(r3)) and showed that it does not have 
a purely fast decay. To understand this point, consider 
the equal time correlation function {Qa{t)Qfiit)Q-y{t)) = 
{QaQpQj) ■ Inserting Eq. (^sj), taking the hmit t oo, 
and noting that in that limit G(t) —^ 0, one obtains 

{QaQpQi) = lim / / / Gar,{t- Ti)Gj3s{t- T2) 

*^°°Jo Jo Jo 

xG^cit " T-3)((/)^(ti)(/)5(t2)(/)c(t3)) 

dTidT2dT3. (40) 

If '/>(t) were Gaussian with zero mean, the three point 
correlation on the right hand side would be zero, but 
since the left hand side of Eq. (40) doesn't vanish, (/)(t) 
is not a Gaussian fluctuating force. 

The three point correlation function in Eq. ( |40| ) can't 
have a purely fast decay either, since, by isotropy, 
((/)^(Ti)(/)5(T2)(/)^(r3)) = 0{k'^), whereas the left hand side 
is 0{k^). It therefore follows that upon integration of the 
slow part of {(j3rj(Ti)(j3s(T2)4>c,{T^)) , one needs to generate 
a factor fc~^. Schramm and Oppenheim obtained the 
explicit form of the slow behavior of the three time cor- 
relation function of the fluctuating force for the case of a 
single slow variable A. It was shown that for r2, T3 > ri 
and small |ki| (other cases are similar), 

('/'ki(Tl)0k2(T-2)0k3(T-3)) 

« 2|kipi^e-^l''^l'l^^-^^l(ik,0k.(r2 - r3)<^k3) 

where D is the diffusion constant, which is clearly a 
slowly-decaying function of T2 — ti. On the other hand, 
it appears that there is a fast decay in T2 — T3, so that 
the whole expression is small when that time difference 
becomes large. Nonetheless, the slow behavior in T2 — ti 
will bring about a factor of 0{l/k^) upon integration. 
For the full restoration of the 0(fc°) term on the left 
hand side of Eq. ( ^ ) , we refer the reader to the original 
paper p5[ . 

Apparently, the assumption that V projects out all the 
slow behavior is not sufficient to specify when a correla- 
tion function is fast-decaying since clearly (j)a (t) itself is 
not a fast variable in every context. Notice that despite 
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the appearance of slow behavior in the correlation func- 
tion ((/)ki (Ti)(/)k2 {T2)(t>'k3 {t3))j there are instances in which 
a multi-time correlation function of fluctuating forces is 
certainly small; namely, when at least two of the time ar- 
guments of the forces are well-separated. If the time ar- 
guments of the fluctuating forces are not well-separated, 
slow behavior can occur. In Eq. (p6|), we cannot assume 
that the integrand is peaked around (r2,ri) = (0,0) be- 
cause these two variables can come arbitrarily close in 
the integration over T2 and ti along the line ti « T2 and 
slow behavior can be expected. We conclude that the lo- 
cal time dependence of the instantaneous Eq. (|3^) cannot 
be justified. 

We propose the following general rule to determine 
when a correlation function is fast decaying: in a corre- 
lation function involving fluctuating forces, the function 
decays quickly in a pair of time arguments, provided these 
are well-separated in time. Note that "well-separated" 
here means that the time difference is larger than the 
microscopic time r^. In applying this rule to situations 
when integrations are carried out over the time argu- 
ments of fluctuating forces, we require that the time ar- 
guments can only get close at isolated points which give 
contributions of measure zero to the integral. 

With these rules in mind, consider the correlation func- 
tion 

(0q(-T-2)0/307(ti)). 

This is clearly fast in T2 as well as in ri provided T2 and 
Ti are positive. Therefore, Eq. should be correct, as 
the correlations of the fluctuating forces in Eq. (^) are 
of the form above, which we can write as 



We observe that the correlation function above is fast in 
Ti and T2 provided ri ^ and T2 3> since it has 
the form of a succession of two fast forward propagations 
which yields an expression in which all time arguments 
are well-separated. If one of the times were negative, i.e., 
one of the propagators propagated backward in time, the 
expression would no longer be (purely) fast. Hence an 
alternative way of identifying terms which are fast in all 
time arguments is to require that they have only forward 
fast propagation when applied in succession. 



B. Correlation functions involving multiple times 

Using the conclusions of the previous section, we will 
now derive a recursion relation for multi-time correla- 
tion functions. We consider the general case, denoted by 



) or G^"''^j({ti}), defined as 



G 



(n) 
{«.} 



with all with i = \ . . .n and ti > so the arguments are 
time-ordered. We write 

(Q5o<3"-> {ti+t2 + ---+ tn) ...Qa2 {h + t2)Qa, (tl)) 



Qa2 



B) 



where B — Qa,^-iQa„{tn) and X is the operator 

X — Q ao^ Qai^ Qa2 ' ' ' ^ 

Using Eq. (^) applied to B, and inserting the result 
into (Xe^*"-iB), we get 



= (g^B)i^;/(XQ^-(t„_i)) 

(5(5 







(Ql e^^^i7'iB>A7/(XQ^(t„_i - ri))dri. 



(Xe^-J-*"-i7'±B) (41) 



From the discussion in Sec. [II A, it is now apparent that 



the third term can be considered fast in because 
only forward propagation occurs in X since all ti are 
positive, and there is a projected propagation in 
For macroscopic times for which tn-i S> t^, this term 
can be neglected. Note, however, this term could not be 
be neglected in integrals of Eq. (^T|) over the time i„_i. 
Inserting B = C - • - - - — 
obtains 



-iQa„{tn) and using Eq. (M), one 



G(n— 1) /, , 
(5,a„_2...v'-"-l ^ '''litn-2, 

-\- fast term in 

where 



drdri 
(42) 



(5(5 



(43) 



Eq. ( P4 ) is the desired recursion relation. Neglecting 
the term which is fast in t„_i, the n multi-time corre- 
lation functions can be related to the n — 1 multi-time 
correlation, and thus ultimately in terms of Ga^aoiti) — 
Gaiaoiii)- The instantaneous version of Eq. ( p2[ ) is sim- 
ply 

^i"^ {tn, ■ ■ ■) ~ Ga„p{tn)Mpa„-iSG'f'^ {tn-1, ■ ■ ■), 

(44) 

where = dri dTMsw{T,Ti). 

Applying Eq. (^ ) to the three time correlation func- 
tion Ga'Yi3{t2,ti), Eq. (^) is recovered: 



Gajp{t2,ti) = Gas{t2)Ms^eGgp{ti). 



(45) 
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Again using the recursion relation, one can also derive 
equations for correlation functions involving four or more 
times, e.g., 

'-^0/37.5(^3,^2,^1) = GaQ{h)Mc^(}gGg,^{t2)Mrj-f\G\s{ti). 

(46) 

It is clear now that any multi-time correlation G^"^ can 
be written as a product of n factors of Gn and n — 1 
vertices M . 



C. TV-ordering of triple time correlation function 

In this section, we consider the leading iV-order and 
mode coupling terms expressions for the triple time cor- 
relation function for linear densities in which \a\ — \ j3\ = 
I7I = 1. Taking Eq. ( ^ ) for the linear densities, this 
correlation function is given by 



^01/3(^2,^1) — Gas{t2)Msi0Ggp{ti). 



(47) 



From this relation, it is evident that the A^-ordering of 
M^i5 follows from that of G~^(i2)Gai/3(t2, ti)G^/(ti), 
and hence we need to establish the iV-ordering properties 
of {Qa(t2 + ti)Qi{ti)Qf3) . One can show, by induction 
in that 

r 0(TVl"l+i) if |a| < |/3| 
{Qc.{t2+h)Qi{h)Qp)^l 0{N\P\) if|a|-|/3| . 

[ 0(Arl/3|+i) if > 

Combining with Eq. (||) and Eq. (|^), we obtain 

r C)(Ari-(l/3|-l"l)) if ^ 

Goi/3(t2,ii) = < 0(A^") if \a\ = |/3| . (48) 

[0{N^) if|a|>|/3| 

Using Eq. ( pl| ) and Eq. (^) to establish the A^-order 
of G;^i(t2)Gc,i/3(t2,ti)G;^5(ti), one finds that the A^- 
ordering of M^is follows the same A^-ordering rules as 
G^is{t2,ti)- 

With the A^-ordering expressions above, the dominant 
contributions to Giii(i2,ii) are given by 

GiiiihM) =Gii{h) * Afin * Gii(ii) 

+Gi2(t2)*M2ii*Gii(ii) (49) 
+Gii(t2) * Mn2 * G2i(ii) + 0{N-^). 



J 



The leading mode-coupling corrections to Eq. ( ^ ) in- 
volve a large number of terms of order A^^^. Collecting 
these additional terms gives a net factor of M w Vk^, 
so that the sum of all A^~^ terms gives a term of or- 
der M /N « (fccC)^ which survives in the thermodynamic 
limit. To first order in M/N, one obtains the following 
correction terms: 

Gi3(i2) * M312 * G2i(ii) + Gi3(t2) * A/311 * Gii(ii) 
+ Gi2(t2) * Af2i2 * G2i(ii) + Gi2(t2) * M213 * G^iih) 

-|-Gii(t2)*Mii3*G3l(tl) 



These results can be easily extended to higher order 
correlations G^^^.-j with n > 3 and where \ai\ — 1 since 
all the necessary A^ orderings are known. Thus, any 
multi-time correlation function of Qi's can be expressed 
in terms of the vertices Maip and two-time (but possibly 
multi-point) correlation functions. In turn, the multi- 
point correlation functions can be expressed in terms of 
the vertices Ma/3 and the linear propagator Gii(i) as 
explained in Sec. [ID, and hence all time dependences 



in multi-time correlation functions at long times can be 
expressed in terms of the two-time correlation functions 
of linear densities. From the mode-coupling formalism, 
the two-time correlation functions can be evaluated self- 
consistently via the relations Eq. (pj) and Eq. ( p3| ) . 

Using the results from Sec. p^I D| , and Eq. (Efl) and 
Eq. ( pO[ ) , one can obtain the leading A^-order expressions 
for Gill in terms of Mu, Mm and Gii(i) as follows. In- 
serting the reduced forms of the vertices M211 and Mi 12 
that are derived in Appendix ^ into Eq. (|4^) yields. 



G'k-q,q,k(^2,tl) — Hi 



H2 + , 



(50) 



where 



Hi - Gk^q(t2) ■ A/k"q;q:k ' ^^,^(^1) 



GlHh 



K. 
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H2 - Gk_q:_q^k(*2) 
Hs = G^_^{t2) ■ Gk^_q q.k(tl) 



In i/2, we can use Eq. ( ^ to express G^_^._^ ]^{t2) in 
terms of the linear Gii(i) and vertices M. For G21 in H-^, 
using the fact that Go,f}{t) = G p^{t)Ks,^Kj^^ [Eq. (H)], 

an analogous expression for G21 is obtained. 



G^'-q,q;k(^l) = S K^' ■ Gi'in) ■ Mlfk-q^q : { K^tl - Tl ) • Xq] O [Gil^{tl - Tl) ■ ATk-q] } dn 

Thus, in Eq. (pO|), we have 

H2 = fj GI\{t2) ■ M^^q;_q,k : { [G"q(t2 " T2) ■ A'q] 

H, = S J^' K^' ■ GI\ti) ■ M^fk-q,q : { K'ih - n) 



o [Gl'{t2-r2)-Gl'{ti)]]dT2, (51) 

^q] ° [Gl\-i.,iti - Tl) ■ Gi\{t2) ■ i^k-q] ]dTi. (52) 
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In a subsequent publication pa^ , we demonstrate that the expressions above give excellent results for a moderately 
dense hard-sphere system in the hydrodynamic regime. 



IV. COMPARISON WITH KAWASAKI'S 
THEORY: NON-GAUSSIAN EFFECTS 

Roughly 20 years ago, Ronis examined higher- 
order correlation functions within the Kawasaki mode- 
coupling formalism ps] , ^ . The treatment itself is too 
technical to recapture here, so we will simply state the 
results from that paper to compare with those from the 
present theory. For a multi-point correlation, Ronis ob- 
tains [his Eq. (3.11)] 

Jo 

xG/3iQi(ki,ti)G'/32Q2(k2,ii) (53) 

We keep Ronis's notation here, as it is close enough to 
ours to be understood. Eq. (|5^) looks essentially like 
Eq. (^0|), but is missing the first term in Eq. ( |3l| ) as 
might be expected from a Gaussian theory for three-point 
correlation functions. Intriguingly, this first term is all 
one would obtain for the three-point correlation function 
from a projection operator approach if only linear densi- 
ties were included in the basis set for the long-time dy- 
namics. In addition, the vertex V in the Gaussian theory 
also differs from the vertex M21 due to the subtraction 
terms in the basis set which are not present. These differ- 
ences in functional form of the vertices can significantly 
alter the time profile of both multi-point and multi-time 
correlation functions 

In Ref. the following expression for a three-time 
correlation function was derived [Eq. (6.5b) therein], 

{to)Al^^,^ (0))) {Ak, AIJ-^^^ 



X G,32a^ (^2 , Ti ) {Ai3^ (n )^ai ,fci (t 1 )) 



+2 



X Gf32a^{k2,Ti){Ap^^^ko{n)Aa„,k„{to))- (54) 

Comparison between Eq. (|5j) and the expression for 
the three-time correlation function in Eq. (|5^) is facili- 
tated by noting that in the instantaneous approximation 
and neglecting mode-coupling corrections, we can write 

Gi\t2 - T2) ■ G]}{h) = G]}{h +t2- T2); 

G'k-q(i2) ■ GiL^ih - Tl) = GkLq(<l + t2 - Tl) . 



From careful inspection of Eq. (|5J) and Eq. { pO\ j , one sees 
that H3 (see Eq. (^)) is essentially equivalent to the sec- 
ond term Eq. (M) . However the term H2 differs from the 



r 



first term in Eq. (54) in two ways. First, the way in which 
the indices are contracted with the vertex Vf°^^J^^i. , as 
written in Ref. [Q, differs from the tensor contractions 
in H2- Second, and more importantly, there seem to be 
significant differences in the upper limits in the time- 
convolution integrals, which in Ref. [Q is to = ti +t2, as 
opposed to t2. This is particularly intriguing in light of 
the observation that the upper limit of to was obtained 
in Eq. ( |37| ) where the time dependence of the correlation 
function of the fluctuating forces was treated incorrectly. 
Nonetheless, in both mode-coupling theories, the higher 
order time correlation functions are expressed in terms of 
ordinary (two) time correlation functions. The major dif- 
ferences between the theories arise because of the Gaus- 
sian approximation in the Kawasaki formalism. Hence it 
is not surprising that some (static) three point correla- 
tions are missed since they vanish if the linear densities 
Qi are assumed to obey Gaussian statistics at all times. 
This deficiency was noted by Ronis who suggested that 
these differences result in significant deviations only at 
short times. It is clear from the present formalism, how- 
ever, that this is not the case since terms of the form 
^11(^2) * Mm * Gii(ti) decay slowly in both ti and t2. 
These findings have been confirmed in numerical simula- 
tions of hard sphere systems [m] . 



V. HIGHER-ORDER CORRELATION 
FUNCTIONS FOR THE IDEAL GAS 

In this section, the mode coupling formalism will be 
illustrated for an ideal gas system in the grand canonical 
ensemble composed of particles of mass m in a volume V 
at an inverse temperature f3. In the ideal gas, the motion 
of each particle j is given by 



r,(t)=r,(0) + 



Pjit) =Pj(0). 



Given the simple form of the particles trajectories, any 
time correlation function can be calculated exactly and 
compared with the expressions that follow from mode 
coupling theory. 



A. Conserved quantities 

An essential step in applying the formalism to a par- 
ticular system is the identification of the slow variables of 
the system. In any gas composed of point particles, par- 
ticle number, momentum and energy are conserved and 
hence their corresponding densities, are slowly-varying 
quantities. The ideal gas system is quite different from 
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simple gases in that it has many more conserved quanti- 
ties since the momentum pj of each particle is conserved 
along all directions. Consequently, a tagged particle den- 
sity of the form pYx'pTv P]z^ ^^^^^ ' where rrix, my and 
are arbitrary integers, should be included for each parti- 
cle j in the set of slow variables. However, for collective 
modes, it is not hard to show that it suffices to include 
densities of all analytical functions of the momenta f{pj ) , 
i.e., 



N 



since the contribution of a single tagged particle to cor- 
relation of extensive variables is 0{1/N). 

Taking the Hermite polynomials as a basis for the 
functions f{p), where 



the complete set of linear slow variables is given by 



3 



(55) 



where {i\ denotes the set of three indices {ixiiy,iz\, 
each of which runs from zero to infinity. If we define 
= p|V/3/2m, = and u| = p|V/3/m, 

then the inner product of the Hermite polynomials cor- 
responds to the canonical average. 



Since (H^u)) — unless i — 0, A^'^ is given by 



(56) 



J\k\pt/m^-Pp^/i2m) 

v/27rm//3 



-dp 



g-|kPtV2™/3 ^ (57) 



where k = \]s.\/y/mP is a conveniently scaled wave- vector. 

We will compare the exact result Eq. (^^ to the re- 
sult from the mode-couplin g fr amework of this paper ob- 
tained using Eqs. (||) and (|lO|). The first point to note 
is that Qa{t) is proportional to Q^, which follows from 
the facts the all Hermite polynomials have been included 
in the set of slow variables and that {d/dt)Hn{u) = for 
all n. This, in turn, implies that V±Qa = 0, and hence 
the fluctuating force (pg (t) vanishes for all t. Accord- 
ing to Eqs. (|l0| ) and (M ), Mapir) = 2S{T)Ma0, since 
4'ait) = 0- Thus Eq. ([13|) is exact, with 



M, 



(58) 



Since Qi is can be written as a linear combination of 
Qi's, and, in general, dotQa = Y,\p\<\a\ "'apQp, we con- 
clude that {QaQ*p) = for (3 > a since the multi-linear 
basis set is orthogonal in mode-order by construction. 
Similarly, since {QaQp) = -{QaQ*p), (QaQp) = for 
a > (3 and hence Map is diagonal in mode order, im- 
plying that Eq. ( [Tsl ) decouples at each mode order n 
into equations for the multi-point correlation functions 

Gnn{t)- 

Focusing on the linear variables, Eq. (|l^) reduces to 



G 



{i}{nQ{i}{j} 



it) 



{1} 



Due to the orthogonality of the ^k^'s, G^*^"^(0) = 
S{i}{j}, and hence. 



G«^^">(i)=[exp(M^H) 



(59) 



and the correlation function {A^^ A^^ ) is given by 



(i« A^^ ) = {N)S,^,A^,J, 



The structure of the matrix Mn is relatively simple, 
which makes it possible to actually calculate Gii(t) from 
Eq. (|59|). Fixing the direction of k to be along x, for 
|k| 7^ we note that. 



B. Two-time, two-point correlation functions 



Ai'^ ^ i~k 



Jx 



Consider the density mode Nf^, which corresponds to 
^k°^ ({0} is short for {0,0,0}). For k 7^ 0, the density- 
density time correlation function is given by 



/^N ,k-(r,(t)-r,)\ 

t-k W - " 



This result follows from the fact that (X^t^i exp(ri -k)) = 
{N)Sko, and the statistical independence of particles. As 
the momentum pi is Gaussian-distributed, one obtains 



where Eqs. ( jSq ) and the recursion relation for the Her- 
mite polynomials 

2uHn-i - 2{n - l)Hn-2 = Hn{u), 



have been used. With Eq. (jSg), this leads to the tri- 
diagonal form for the matrix Mn at linear order. 



M, 



{i}{3} 



:{\/ix + 1(5^^+1 + \fQ5i^-ij^). (60) 
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The terms in Eq. ( |6C| ) can be conveniently expressed as 
the expectation values of the raising and lowering opera- 
tors in the Hermite basis representation of the quantum 
Harmonic oscillator. Identifying 



The time correlations for other basis functions are now 



= A 



{n,jy,j^} 



and defining 



{m\Bn) = {Ai"'''-'''^* BAi^-'--''^)/{N) 
for any operator B, and 

a\n) = \fn\n — 1); a}\n) = Vn + l\n + 1), 
one has 

(77l|a7l) = 



In this representation, the matrix Mj^ 
as 



{m\a^ n) = Vn + 1 Sm,n+i- 

^'^^■'^ can be written 



where A4 is the operator 

M = ik{a^ + a). 



(61) 



In essence, A4 is nothing but the Liouville operator C, 
restricted to act on the space of phase space functions 
which are linear combinations of Ak- 

To reproduce Eq. (|57|), according to Eq. (p9|), we must 
evaluate 

(0|e-^*0), 

which can be done in straightforward fashion using the 
Campbell-Bakcr-Hausdorff formula [Q. This formula 
states that if A and B are linear operators, a linear op- 
erator C exists such that e^e'^ = e*", where C of the form 
C = A + B + ^[A,B] plus repeated commutators. Tak- 
ing A = ikta^ and B — ikta and noting that [a, a^] = 1, 
Eq. (|l|) yields. 



which can be re-arranged to give, 



(62) 



and therefore 



(0|e-^*0) = (e-*'''*''0|e**=*°0) 



^{ktf/2 



But as e*'=*"|0) = |0), 

(0|e^*0) -e-(^*)'/^ 

which coincides, as expected, with the exact result 
Eq. (§). 



easily obtained using Eq. (62). For example, 



(N) 



(63) 



Note that it is only with an infinite set of conserved 
quantities that it is possible for a system with no dis- 
sipative terms to lead to relaxation. One might argue 
that taking a finite number of slow variables, (say, mass, 
momentum and energy density) the dissipative terms in 
M would no longer be zero. While this is correct and 
M is not purely imaginary, M would still be finite in 
each mode block, and one would get exponential decay 
in the instantaneous limit instead of Gaussian decay for 
all correlation functions. 



C. Two-time, multi-point correlation functions 

Consider the correlation function 



^{0}{0}{0}^^^^ ^ 
^k-q,q;k 



with |k| 7^ 0, |q| 7^ and k 7^ q. The direct calculation 
of this correlation function is simple if one notes that 



N 



N 



1=1 



(k-q)-n(t: 



N 
n=l 



■•")/(iV) 



N 

1=1 



since any term in the summation over particle indices in 
which I ^ m, I ^ n OT m ^ n yields a Kronecker delta for 
one of the wave- vectors. Such terms do not give a contri- 
bution if all wave-vectors are non-zero. The wave-vector 
q drops out of the expression, and 

N 

^.iT<J;rW = (E^*-"'^'")/W 

-{ktf/2 



1 = 1 
i{0 



(64) 



The mode coupling derivation of this same result goes 
as follows. From Eq. (pQ), noticing that M21 — 0, we 
immediately obtain G2i{t) = 0. The simplest argument 
for the vanishing of G2i{t) is that since M is diagonal 
in mode order, no correlations between different mode 
orders exist, so that 

G^p{t) = ii\a\^\(3\. 

From the definition of the multi-linear basis set in 
Eq. (^) , cj^_}^^^.^^ (t) can be written in terms of the lin- 
ear and bi-linear densities as 
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However since G2i{t) — 0, C^'i'tj q^k°^(*) S^*^ i^s value 
solely from the subtraction terms. Furthermore, since 
{Nk-qNqN^) = {N) for an ideal gas system, one gets 

in agreement with Eq. (|6^). 



D. Multi-time correlation functions 

Finally, we conclude our discussion of the higher-order 
correlation functions in the ideal gas by examining the 
validity of the expressions for the multi-time correlation 
functions in Eq. ( [45| ) . For the three time correlation func- 
tion for the density mode Nk, direct calculation gives. 



^k-q,q,k (t2,tl) - { 



N 

g»Pr[k(ti+t2)-qt2]^ 

1 = 1 



^ g-|k(ti+t2)-qi2r/2 



(65) 



where k = k/ ^/mj3 and q = q/ ^/m(3. 

If k and q are both along i;, the mode-coupling expres- 
sion for the multi-time correlation [see Eq. (Eq)] is 



G 



{0}{0}{0} 



(i2)M, 



k-q,q,k ^k l^^lJ' 



where the repeated sets of indices {«} and {j} are 
summed. Since G^'^^^{t) was previously evaluated in 
Eq. (H), and G{^>^(t) = G^^'J* {-t), the only unknown 
quantity in this expression is 

Writing this out using Eq. ( p5| ) , the only surviving terms 
in the summations over particle index for an ideal gas 
system are the terms where all indices are equal, so 

^"k-q,q,k —"b}b}- 

Hence 

oo 



1 



J. 



This can be summed to 

p-[|k-q|"t^-|k|"t?-fc4L-gx)tlt2]/2 



which corresponds to the exact result Eq. (pq 



VI. SUMMARY 



In this paper, a mode coupling theory was presented 
in which multi-point and multi-time correlation functions 
are expressed in terms of ordinary two-point, two-time 
correlation functions and a set of vertices. The mode cou- 
pling theory developed here does not assume that fluctu- 
ating forces (noise) are Gaussian distributed and there- 
fore is a generalization of mode coupling theories based 
on Kawasaki's formalism Furthermore, unlike 

kinetic theories, it is not restricted to low densities and 
should be applicable to dense fluids where cooperative 
motions of particles and collective modes are important. 

The formalism is based on projection operator tech- 
niques, which, for ordinary two-point, two-time correla- 
tion functions, lead to a generalized Langevin equation in 
which the memory function decays on a microscopic time 
scale. The simple extension of the projection operator 
formalism to multi-time correlation functions is compli- 
cated by the fact that the fluctuating forces appearing in 
the generalized Langevin equation do not obey Gaussian 
statistics. Furthermore, multiple time correlations of the 
fluctuating force can in fact have a slow decay when the 
time arguments of these forces become comparable. 

In order to treat multiple time correlation functions 
of fluctuating forces properly, the correlation functions 
were massaged so that the time arguments of all fluc- 
tuating forces appearing in the correlations were guar- 
anteed to be well-separated, ensuring that all memory 
functions which arise in the mode coupling theory de- 
cay to zero on a molecular time scale. This construc- 
tion allows equations which are local in time to be ob- 
tained which relate the multi-time correlation function to 
two-time but multi-point correlations coupled by essen- 
tially time-independent vertices. The multi-point corre- 
lations, in turn, can be written as convolutions of two- 
point and two-time correlation functions coupled by time- 
independent vertices. These correlation functions can ei- 
ther be taken directly from experiment, simulation, or 
can be solved self-consistently within the mode coupling 
formalism. The vertices, which are composed of a static 
part (Euler term) and a generalized transport coefficient, 
can similarly be calculated from kinetic theory or taken 
from molecular dynamics and Monte-Carlo simulations. 

The equations for higher-order correlation functions 
contain an infinite sum of terms which can be made 
tractable for systems with a finite correlation length by 
applying a cumulant expansion technique, termed the A^- 
ordering method. The method was applied to obtain the 
leading order and first mode coupling corrections of ex- 
pressions for multi-point and multi-time correlation func- 
tions. A key step in the A^-ordering method and the 
proper setup of the theory is the definition of an orthogo- 
nal multi-linear basis. Although in principle it is possible 
to apply cumulant expansion methods to other choices of 
a multi-linear basis, the orthogonalization procedure sim- 
plifies the perturbation analysis enormously and helps to 
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avoid erroneous truncations of the mode coupling series. 

The expressions for the higher-order correlation func- 
tions bear a resemblance to those found by Ronis |^] 
within the framework of Kawasaki's mode coupling the- 
ory. In this approach, the linear densities composing the 
set of slow variables are assumed to be Gaussian random 
variables at all times. Although the Gaussian assump- 
tion provides another method of simplifying the mode 
coupling series for higher-order correlation functions, cer- 
tain terms are absent from the Gaussian theory which are 
neither small nor quickly-decaying. 

The mode coupling predictions for higher-order corre- 
lation functions for an ideal gas system were calculated 
analytically and shown to give the exact results for both 
multi-time and multi-point correlation functions. An es- 
sential step in arriving at the correct result was the in- 
clusion of a complete set of densities in the set of slow 
variables. Although the ideal gas system does not con- 
stitute a rigorous test of the formalism since all fluctuat- 
ing forces vanish, it is important to note that the formal 
mode coupling theory expressions for the higher-order 
correlation functions yield the exact result. In a future 
publication p^ , we compare the mode coupling predic- 
tions for the multi-point and multi-time correlation func- 
tions for a hard sphere fluid to data from simulations. 
The theoretical predictions of all higher-order correlation 
functions are in remarkable agreement with the simula- 
tion results in the hydrodynamic regime provided both 
Eulcr and dissipative vertex couplings are included. 

The theory outlined here has obvious applications to 
multi-dimensional Raman and NMR spectroscopy, and 
simulation studies of dynamic heterogeneity in dense flu- 
ids, glasses and polymers. Since the theory involves phys- 
ical correlation functions, it is well-positioned to address 
fundamental issues in characterizing the dynamics in sys- 
tems exhibiting non-exponential relaxation processes and 
frustration. In fact the current formalism has been used 
in Refs. [E6 52| to justify some of the approximations 



made in mode-coupling theory for the supercooled liquids 
pq,Bl|. These avenues are currently being investigated. 
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APPENDIX A: REDUCTION OF THE VERTICES 

M211 AND M112 

Neglecting mode coupling corrections, the vertices 
M211 and M112 can be reduced to very simple forms. 
The strategy used to simplify these terms is similar to 



that used for the factorization of Mq^ in Ref . Q| . First 
M112 is re-written with the help of Eq. ( ^5|) as. 

Using the A^-ordering of Gaii3{t2,ti), 6*10(^2), and 
G/32(ii), one sees that to leading A^-order, 



M1I2 = G^^'{t2)Gii2it2,h)G22\h). 



(Al) 



In this expression, the leading A^-order contributions are 
obtained from the part of the various factors in Eq. ( Al) 
which are diagonal in wave- vector. 

By the property in Eq. ([l^), Ar22 can be factored as, 

{At_^iti+t2)A'^{h)Qt^,^:) 

In a similar fashion, K22 factors to leading A-order, 



-"^k-q.qik-q'.q' ~ ^k-q " ^q "qq 



k-q- 



i^-'^k- 



k-q,q', 



where the notation introduced in Eq. ( p9| ) has been used. 
This leads to the following factorization of the wave- 
vector-diagonal part of Gii2{t2,ti): 

G'k-q;q;k-q',q'(*2,tl) ~ G^-q{h + ^2) ° (^1 )<5qq' 

+ Gj,Lq(il+i2)»Gll(il)'^kq-q'. 

The factorization of G^} can be worked out as well: 



q,q;k— q',q' 



+ Gk-q (tl) • (i^f^kq-q' 



Inserting these expression into Eq. (Al) and using 
Eq. (|T7|), yields 



k— q;q;k— q',q' 
« G^Lq (t2) • GlLq(t2 +h)- GllTq {h) O Ijqq, 

+ Gl':^ it2) ■ Gll_q(i2 + h) ■ Glljq ih) . l<5qq, 

W lol,5qq, + l.l,5k-q,q', (A2) 

where l""^ — Sac- 

Finally, Eq. (A2) and the relation M211 = Ar22 * 
{Mil* 2)* * ^11 gi'^^ the leading A^-order expression for 
M211, 



k-q'+q,q';-q;k ~ ^q O 1 (^kq' + 1 • A'„ 5, 



q "qq' > 



(A3) 



to lowest order in the mode-coupling parameter M/N. 
It interesting that the first term in Eq. ( ^ ) gives the 
only contribution to M211 at order {M/Ny: In other 
words, when mode-coupling corrections are neglected. 
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the term involving the fluctuating forces in Eq. ( p3| ) can 
be dropped. 

Note that in this appendix, we have allowed ourselves 
to neglect the difference between the inverse of quantities 
on the multi-linear level and the linear-linear sub-block 
level, which is correct to order {AI/N)'^ and consistent 
with the lev el o f approximation of the rest of the deriva- 
tion of Eq. @ and Eq. @. 



APPENDIX B: SYMMETRY PROPERTIES OF 

Gc,f,{t) AND M^f,. 

In the ideal gas case, the Qa's are either even or odd 
functions of momentum. This property implies that the 
elements of the multi-linear basis set Qa are either sym- 
metric or anti-symmetric under the time-reversal oper- 
ator T, which reverses the momenta. Mathematically, 
this property can be written as TQa — jaQai where ja 
is either 1 if a contains an even number of momenta in- 
dices, or —1 when it has an odd number. In addition, 
there typically is a symmetry under the reflection opera- 
tor TZ which inverts both the momenta and the positions 
of all particles. As the basis set elements Qa depends on 
the spatial degrees of freedom through exp(jk • r^) [see 
Eq. (|55|)], it follows that TZQa = laQ*a- These relations 
also holds for other systems in which the potential en- 
ergy depends only on the distances between particles in 
the system. These two symmetries plus time translation 
invariance have the following implications for Map and 
Gap{t) (see also Ref. H). 

As the equilibrium distribution function is invariant 
under TZ, 

Gap{t) = {Qa{t)Ql)K^^ = ((7^Q„(^))(7^Qp>i<^^^ 

^IclpGlpit). (Bl) 

Hence, if 7a 7/3 = 1, the imaginary part is zero, otherwise 
the real part is zero. Since the wave-vector dependence 
of the densities always enters in the form of i times a 
wavevector, imaginary correlation functions must be odd 
functions of the wave- vector and real correlation must be 
even functions of wave-vector, provided these quantities 
are analytic in wave- vector. 



Time reversal invariance, Te 



ct 



-a 



T, yields 



((e^*Qo)Q;) = {{Te-^'TQa)Qp) = {{e-^'TQa)TQ*p) 

and hence Gagjt) = ja^pGapi—t). This can be com- 
bined with Eq. (Bl) to yield 



Gafl{t) = Glp{-t). 



(B2) 



So from Eqs. ( p31| ) and (p32D, we conclude that if 
7a7/3 = 1) Ga/sit) is real, even in wavevectors and sym- 
metric under t —>■ —t, whereas if 7a 7/3 = —1, it is imagi- 
nary, odd in wavevectors and anti-symmetric under time 



reversal. As M = {dG{t)/dt)G~^{t), this also implies 
that Ma/3 is real and even in wavevectors if 7a 7/3 = 1, 
and imaginary and odd in wavevector otherwise. 

The following ordering is now valid when the magni- 
tudes of the wave-vectors are small. Imaginary correla- 
tion functions and vertices, being odd in the wavevec- 
tor arguments, are typically of linear order in the wave- 
vectors. But as the vertices M contains time derivatives, 
they will always be at least of the order of the wave- 
vectors, so a real-valued vertex is at least quadratic in 
wave-vector, whereas a real-valued correlation function 
is typically of order one. 

Finally, using time-translation invariance, Gap{t) = 
G*~(-t)KaaK~l, which conveniently combines with 
Eqs. (^) to 



Gafi{t)^Gp^{t)KaaK-l. 



(B3) 
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